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Vibrationally excited CO2, formed by two-body recombination from CO(*X*) and 
O(?P) in the wake behind spacecraft entering the Martian atmosphere, is believed to 
be responsible for the higher than anticipated radiative heating of the backshell, 
compared to pre-flight predictions. This process involves a spin-forbidden transition 
of the transient triplet COz molecule to the longer-lived singlet. To accurately predict 
the singlet-triplet transition probability and estimate the thermal rate coefficient of 
the recombination reaction, ab initio methods were used to compute the first singlet 
and three lowest-energy triplet COz2 potential energy surfaces and the spin-orbit 
coupling matrix elements between these states. Analytical fits to these four potential 
energy surfaces were generated for surface hopping trajectory calculations, using 
Tully’s fewest switches surface hopping algorithm. Preliminary results for the 
trajectory calculations are presented. The calculated probability of a CO(‘Z*) + O(°P) 
collision leading to singlet COz formation is on the order of 10-4. The predicted 
flowfield conditions for various Mars entry scenarios predict temperatures in the 
range of 1000K-4000K and pressures in the range of 300-2500 Pa at the shoulder 
and in the wake, which is consistent with a heavy-particle collision frequency of 10° 
to 10’ s+. Owing to this low collision frequency, it is likely that CO2(42,*+) molecules 
formed by this mechanism will mostly be frozen in a highly nonequilibrium ro- 
vibrational energy state until they relax by photoemission. 
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L Introduction 


In order to simulate the flowfield around a spacecraft and design efficient thermal protection 
system for entry into a planetary atmosphere, one needs to understand non-equilibrium chemical 
kinetics and radiation phenomena in hypersonic flows. Work is underway within NASA to develop a 
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non-equilibrium model based on ab initio computational chemistry that can be incorporated into 
computational fluid dynamics (CFD) simulations and lead to predictions of the convective and 
radiative heat fluxes experienced by the spacecraft. The chemical and physical processes 
encountered during high-speed Earth entry are fairly well understood, but the situation for Mars 
entry is much less so. 


The Martian atmosphere is composed mostly of CO2 (96% by volume), with some Nz (1.9%) 
and Ar (1.9%) plus other minor constituents. [1] For most proposed NASA missions, the vehicle entry 
speeds will be 4-8 km/s and the equilibrium temperature of the bow shock layer will be 5,000-9,000 
K. For 6-8 km/s entry, most of the COz is rapidly dissociated in the shock layer, resulting in the 
formation of CO(Z*) and O(3P). In this case, CO 4t-positive emission in the VUV is the major 
component in the radiative heat flux. However, for 4-5 km/s entry, hot COz2 is also present in the 
shock layer and infrared emission from that molecule is an important contributor to the radiative 
heating. See Johnston and Brandis [2] for a discussion of the radiative heating. As the flow expands 
around the shoulder and into the backshell region of the spacecraft, the temperature and pressure 
drop, and recombination reactions are more probable. In this part of the flow, temperature will 
probably drop to 1,500-4,000K and pressure will drop from about 50,000 Pa in the bow shock layer 
to less than 3000 Pa. The radiative heat flux incident on the backshell is predominately infrared 
emission from COz and CO. In practice, the Carbon Dioxide Spectroscopic Database (CDSD) [3] is used 
to predict the spectral intensity of this radiation. That model assumes the populations of the ro- 
vibrational levels of COz2 are in thermal equilibrium, which is probably not accurate for expanding 
flows where relative internal energy populations become frozen as the temperature and pressure 
drop. 


In general, the heat flux experienced by the forebody of the vehicle during atmospheric entry 
is much greater than that experienced by the afterbody. Consequently, the forebody heat flux has 
been studied to a much greater extent by flight instrumentation, CFD and ground-based experiments. 
However, a recent CFD study by da Silva and Beck [4] suggested that non-equilibrium COz infrared 
radiation (IR) is the dominant contributor to the afterbody heating of Mars entry vehicles and this 
radiative heat flux is significantly greater than previously assumed. As a result, the requirements for 
backshell thermal protection for Mars entry are being reconsidered. More recent computational [5, 6, 
7] results seem to support this conclusion. Two sets of expansion tube experiments [8, 9] are in 
general agreement with the CFD studies, but their results are not conclusive. These experiments, 
using fast CO2 flows, have provided the first measurements of radiative heating in the expanding flow 
region behind small sphere-cone models. However, the expansion tube measurements were carried 
out for freestream velocities under 5 km/s and not much of the CO; in the flow was dissociated in the 
bow shock layer. Therefore, their relevance to higher speed entries is unclear. 


The most widely used non-equilibrium chemistry model for describing Mars entries is the 
two-temperature (i.e., T-T,,) model of Park et al. [10] (referred to as Park94). That model includes 
collisional dissociation rate coefficients for COz2 and CO taken from shock tube experiments carried 
out between 1964 and 1983. For this model, dissociation rate coefficients are described by an 
average temperature T,,, which is the geometric mean of the translational temperature T and the 
vibrational temperature T,,. In a recent update by Johnston and Brandis [2], some of the original rate 
coefficient parameters from the Park94 model were adjusted to better describe CO and CO2 emission 
spectra measured in the Electric Arc Shock Tube (EAST) facility at NASA Ames. In the older Mars 
chemistry models, little attention was given to CO2 dissociation, because COz was assumed to be 
rapidly dissociated for most Mars entry velocities. In these nonequilibrium chemistry models, rate 


k 
coefficients for reverse reactions are obtained by microscopic reversibility (Ke, = = with k, and k, 
a i 
defined as the forward and reverse rate coefficients, respectively, and K,, is the equilibrium 
constant) so the flow will relax to the correct chemical equilibrium solution. Collisional dissociation 


of COz2 is known to form CO and O in their ground electronic states, 1Z+ and 3P, respectively. Even 
though COz is a singlet molecule in its ground electronic state, it dissociates along a triplet potential 
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energy surface (PES). Therefore, it must undergo a transition from being on a singlet to a triplet PES. 
However, details about this non-radiative singlet-triplet transition are not well known. 


Over the last 50 years, a number of shock tube experiments have been carried out to 
determine the rate coefficient for COz dissociation. Early shock tube experiments were hampered by 
the presence of impurities, especially organic species, which greatly accelerated the dissociation 
process by catalytically lowering the activation energy. [11, 12] The rate coefficient used in the 
Park94 model was based on the shock tube experiments of Davies [13], because that was deemed the 
most accurate for high temperatures (even though it had an anomalously low activation energy 
according to [11]). Park [10] adjusted the parameter C in the Arrhenius expression (Kp(Tgy) = 
-E 
F(M)-C+Tay": e Ray), where kp is the dissociation rate coefficient, f(M) is the relative efficiency of 
the collision partner M in promoting dissociation, Fp is the dissociation energy and R is the gas 
constant) to match Davies’ published rate coefficient at 5000 K while keeping n fixed at 1.5 and using 
the accepted value for the COz2 dissociation energy (526 kJ/mol). If the molecule were to dissociate to 
the lowest energy singlet asymptote, CO(!X*) + O(‘D), Ep would be 715 kJ/mol and the exponential 
term in the Arrhenius expression would by 100 times smaller at T,,, = 5000 K. There have been 
other measurements of the COz2 dissociation rate coefficient published since Davies’ work was 
published. [14, 15, 16, 17] Jaffe reviewed these experiments and discussed the role of low-lying COz 
triplet electronic states in 2011. [18] All of these experiments were carried out for dilute CO2-noble 
gas mixtures, so the observed excitation and dissociation of COz2 was due to collisions between COz, 
and Ar or Kr. For the Park94 model [10] the dissociation rate coefficients for CO2 + M + CO + O used 
f(M) = 1 for noble gases, 15 for C, N and O atoms, and 10 for CO, CO2, N2, O2 and NO molecules. 
When the actual data points from all of these shock tube experiments (not the Arrhenius fits) are 
plotted together [18], the results are seen to be in good agreement, as shown in Figure 1. The linear 
least-squares fit to all the data points in the figure yielded the following expression for the 
dissociation rate coefficient: kp = 3.83 x 10*!4 + exp(-52116/T) cm3molecules:'s-!, where T is in Kelvin. 
The dissociation temperature of 52116 K corresponds to Ep = 433.31 kJ/mol. This rate coefficient is 
approximately a factor of two larger than the corresponding rate coefficient in the Park94 model. 


The objective of the present study is to investigate the fundamental chemical kinetics process of COz 
formation by recombination and determine the radiative flux from the COz formed in the backshell 
region of spacecraft using ab initio quantum chemistry computations and classical mechanics 
simulations. This process has been studied previously at higher pressure conditions where three- 
body collisions stabilize the nascent product molecules. [19] However, at the low pressure conditions 
of interest for the present study, the frequency of three-body collisions will be much lower than the 
frequency of two-body collisions. Therefore, we have chosen to investigate the following 
recombination process: 
CO(12*) + OP) > CO2(2A’, 3A”) > CO2(22°), 

in which CO and O in their ground electronic states recombine to form triplet CO2, followed by a 
transition between triplet and singlet electronic states that ultimately results in the formation of 
singlet COz. In fact, for a typical Mars entry, we estimate the mean time between collisions in the 
afterbody flowfield will be 0.1-0.5 us. Actually, three triplet electronic states originate from the 
reactant asymptote, one of 3A’ symmetry and two of 3A” symmetry. They are referred to as the 3A’, 
13A” and 23A” states herein. Without a subsequent collision with another atom or molecule, the 
nascent triplet COz will have a total energy above its dissociation limit, and will have an extremely 
short lifetime before it dissociates. Two-body recombination cannot occur unless there is some 
mechanism for triplet-to-singlet intersystem crossing. The ground state singlet potential energy 
surface (PES) has a seam of intersection with each of the three triplet PESs. For COz, each PES is a 3- 
dimensional hypersurface (V(17,,72, 0), with r, and r, representing the two C-O bond lengths and 6 
the O-C-O angle), and the intersection seam has two degrees of freedom. Figure 2 shows a cut of the 
four PESs where the O-C-O bond angle is fixed at 110°, one C-O bond distance is fixed at 1.15 A, which 
is close to re, the equilibrium bond length in CO2. The other C-O bond distance is allowed to vary. The 
ground state curve intersects the three excited state curves at three different points. Under the Born- 
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Oppenheimer approximation, the triplet and singlet adiabatic PESs do not interact; thus, the triplet to 
singlet electronic transition is regarded as spin-forbidden. However, the true Hamiltonian contains 
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Figure 1. The measured rate coefficient data points for CO, dissociation from which 
the experimentally determined rate coefficient fits have been determined. Each 
color symbol represents a different experiment. The Park 94 [10] recommended 
rate coefficient is shown as the black dashed line. Experimental data are from 
Davies (stars) [13], Burmeister and Roth (blue squares) [14], Fujii (triangles) [15], 
Oehschlaeger et al. (diamonds) [16] and Saxena et al. (circles) [17]. The meta-fit of 
all the experimental data points [18] is shown as the green dashed line. Rate 
coefficient units are in cm?mol-s-1. 


spin-orbit coupling terms which enable nonadiabatic electronic transitions, also known as 
intersystem crossing, between states of different spin multiplicities. Spin-orbit coupling [20, 21, 22, 
23, 24, 25] is not only important for describing heavy elements, but it is also crucial in light atom 
systems. For example, it enables intersystem crossing and phosphorescence of excited triplet states 
in organic molecules, it alters chemical reaction paths and it leads to the fine structures in high- 
resolution spectroscopy. With spin-orbit coupling, there is a probability for a triplet COz molecule to 
make a transition onto the singlet PES and become a singlet COz. The reverse transition is also 
possible, meaning there is a probability for a singlet COz to transition into a triplet COz. The 
nonadiabatic transition probability is greatest when the spin-orbit coupling between the interacting 
states is strong. A widely-used model for describing the transition probability between two states is 
the Landau-Zener(LZ) model [26, 27, 28] where the transition probability is: 
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( -2mHs9" } 

Piz = 1-—e hlOV3/dq—AV1/dqlv ; 

where Heo is the spin-orbit coupling between the singlet and triplet state, v is the nuclear velocity 
perpendicular to the crossing seam, and |0V;/0q — OV,/0q| is the absolute difference between the 
gradients of the two PESs at the crossing point. The LZ model is based on a highly simplified 
description of the potential energy surfaces. Tully has developed a more general model called the 
“fewest switches surface hopping model” [29] that is better suited for use in classical and 
semiclassical scattering calculations. If the triplet-singlet transition occurs, under the conditions 
being considered, the nascent singlet COz molecule has a ro-vibrational energy that is below its 
dissociation limit, CO(‘2+) + O(‘D), and thus is trapped in a high-lying ro-vibrational level of the 
ground electronic state. It can only relax to lower levels through subsequent collisions or photon 
emission. 


Relative Energy (kcal/mol) 


05 1 1.5 2 2.5 3 


T(0-co) (A) 


Figure 2. Calculated MRCI+Q potential energy curves of CO2 singlet ground state 
and three triplet excited states as a function of OC-O distance while keeping the 
0-C-O bond angle 6 fixed at 110° and the other C-O bond distance fixed at 1.15A. 
The black curve denotes 1A’ state, the red curve denotes 3A’ state, the green 
curve denotes the lower 3A” state and the blue curve denotes the higher 3A” 
state. The singlet curve asymptote is CO(‘X*) + O(‘D) while the triplet curve 
asymptotes are CO(1x*) + O(3P). 


The ultimate goal of the current investigation is to predict the amount of singlet COz 
molecule formation due to the triplet-to-singlet transition, by computing the thermal rate coefficient 
for singlet COz formation, and the distribution of COz internal energy as a function of temperature 
and initial ro-vibrational level of CO. To achieve this goal, the first step is to compute the four PESs of 
COz by solving the electronic Schrédinger equation using highly accurate electronic structure 
methods, locate the molecular geometries where the crossing between 3A’ and 1A’, 13A” and 1A’, and 
23A” and 1A’ take place, obtain spin-orbit coupling matrix elements at the crossing geometries, and 
assess the likelihood of electronic transitions. After obtaining information on the PESs, the second 
step is to run surface-hopping quasi-classical trajectory calculations for CO(x*) + O(3P) collisions on 
the four surfaces. 


These calculated recombination rate coefficients can be combined with flowfields from CFD 
calculations to estimate the radiative heat flux on the backshell of Mars entry vehicles. 
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II. Computational Methods and Results 


A. Calculated 1A’, 3A’, 13A” and 23A” potential energy surfaces and spin-orbit coupling matrix 
elements 


The ab initio calculations of the electronic energy at various geometries were carried out 
using the Complete Active Space Self-Consistent Field (CASSCF) [30, 31, 32, 33] and Multireference 
Configuration Interaction (MRCI) [34, 35, 36] as formulated in the Molpro suite of quantum chemical 
programs. [37] In order to accurately and consistently describe the four surfaces on equal footing, all 
four states are included in the state-averaged full valence CASSCF and MRCI calculations combined 
with correlation consistent basis sets aug-cc-pVTZ [38], which properly account for electron 
correlation and bond breaking. This also enables us to obtain the spin-orbit coupling matrix elements 
between the singlet and triplet states. All calculations are performed with the Molpro Suite of 
quantum chemical programs (version 2010.1). [37] 


We focus on regions of the PESs near the crossing seams and sampled over 9000 points. The 
location and energetics of the crossing points depend on the geometric variables and some of the 
trends we have identified based on our calculations so far are: (1) 3A’ crosses !A’ at smaller CO-O 
distance compared to 13A”, and much smaller compared to 23A” in general; (2) The relative energy 
with respect to the triplet asymptote CO(‘=*) + O(3P) at the crossing point is lowest for 3A’, then 13A” 
and lastly 23A”. This means the trajectories are most likely to have access to the crossing points of the 
3A’ and 13A” surfaces with 1A’ surface. 23A”-1A’ crossing points lie quite a bit higher in energy 
compared to the other two triplet states, and the likelihood of triplet-singlet transition occurring on 
the 23A” surface is thus likely to be very small. In addition, the global minima on the 3A’ and 13A” 
surfaces are bound with respect to the triplet asymptote while global minimum on the 23A” surface is 
not. thus we are going to focus on the 3A’ and 13A” surfaces; (3) As the bond angle @ increases, the 
3A’-1A’ crossing points increase in energy, and the 3A’-!A’ and 13A”-1A’ crossing points become closer 
together. Figure 2 and Figure 3 are two example cuts of the four PESs with the bond angle 6 fixed at 
110° and 140°, respectively. As expected, the 1A’ curve has a much deeper well. Earlier lower-level ab 
initio calculations by Hwang and Mabel [39] reported similar crossing geometries for the lowest- 
energy singlet and triplet electronic states of COz. 
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Figure 3. Calculated MRCI+Q potential energy curves of CO2 singlet ground state 
and three triplet excited states as a function of OC-O distance while keeping the 
O-C-O bond angle 6 fixed at 140° and the other C-O bond distance fixed at 1.15A. 
The black curve denotes 1A’ state, the red curve denotes 3A’ state, the green 
curve denotes the lower 3A” state and the blue curve denotes the higher 3A” 
state. The singlet curve asymptote is CO(1x+) + O(‘D) while the triplet curves 
asymptote is CO(4X*) + O(P). 
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The spin-orbit coupling matrix elements are computed using the full Breit-Pauli operator 
implemented in Molpro [40] with the adiabatic MRCI wavefunctions. There are three Cartesian 
components of the spin-orbit Hamiltonian, LSx, LSy and LSz. For 3A’-1A’ coupling, only the LSz 
component is non-zero, while for 3A”-!A’ and 3A’-3A” coupling, both the LSx and LSy components are 
non-zero. The transition probability depends on the strength of the spin-orbit coupling, and the 
coupling matrix elements depend on the geometric variables. Figures 4 and 5 show the spin-orbit 
coupling matrix elements as a function of CO-O distance while keeping the angle @ fixed at 110° and 
140°, respectively, with the other C-O bond distance kept fixed at 1.15 A, corresponding to Figures 2 
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Figure 4. Spin-orbit coupling matrix elements as a function of OC-O distance while keeping the 
angle @ fixed at 110° and the other C-O bond distance at 1.15 A (the same geometries as in 
Figure 2). The black curve denotes the 3A’-1A’ spin-orbit coupling matrix element LS7, the red 
solid curve denotes the 13A”-1A’ coupling matrix element LSx, the red dashed curve denotes 
the 13A”-1A’ coupling matrix element LSy, the purple solid and dashed curve denote the 23A”- 
1A’ coupling matrix element LSx and LSy, respectively, the green solid and dashed curve denote 
the 13A”-3A’ coupling matrix element LSx and LSy, respectively, and the blue solid and dashed 
curve denote the 23A”-3A’ coupling matrix element LSx and LSy, respectively. 


and 3. The coupling matrix elements vary with the OC-O distance and 3A’ has the strongest coupling 
with !A’ overall. We have focused on evaluating the matrix elements for geometries near the crossing 
seams since what matters the most are the spin-orbit coupling strengths in those regions. Figure 6, 7, 
and 8 plots the matrix elements as a function of the energy difference between the coupling states. A 
small energy range, +1 kcal/mol, is chosen for the plots. For 3A’-!A’, the coupling matrix elements are 
in between 50 and 85 cm". For 13A”-1A’, the LSx component mostly lies between 30 and 40 cm! 
while the LSy component lies in between 0 to 40 cm". The 23A” and !A’ surface cross at about ten 
calculated points, because we are focusing on the crossing between the lower-lying triplet surfaces 
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Figure 5. Spin-orbit coupling matrix elements as a function of OC-O distance while keeping the 
angle @ fixed at 140° and the other C-O bond distance at 1.15 A (the same geometries as in 
Figure 3). The black curve denotes the 3A’-1A’ spin-orbit coupling matrix element LS7, the red 
solid curve denotes the 13A”-1A’ coupling matrix element LSx, the red dashed curve denotes 
the 13A”-1A’ coupling matrix element LSy, the purple solid and dashed curve denote the 23A”- 
1A’ coupling matrix element LSx and LSy, respectively, the green solid and dashed curve denote 
the 13A”-3A’ coupling matrix element LSx and LSy, respectively, and the blue solid and dashed 
curve denote the 23A”-3A’ coupling matrix element LSx and LSy, respectively. 


with the singlet surface, and both the LSx and LSy matrix elements are between 0 and 30 cm. Hwang 
and Mebel [39] computed values for the spin-orbit coupling matrix element of 20.3 cm for the 
lowest energy singlet-triplet crossing geometry in C2y symmetry (OCO bond angle of 105° and CO 
bond lengths of 1.26 A) and 89.5 cm: for a linear geometry with CO bond lengths of 1.13 and 1.90 A. 


B. Analytical representations of 1A’, 3A’, 13A” and 23A” potential energy surfaces 


In order to run trajectories, we need good analytic representations for the four surfaces. For 
the three triplet surfaces, we use the surfaces from our previous work. [41] While those triplet 
surfaces did not include the present data in their determination and the emphasis was on other 
regions of the PES rather than the singlet-triplet crossing seams, they provide a useful base to start 
this study. The predicted shape of the PES in the vicinity of the crossing seams is in good agreement 
with the present calculations, although the overall energy is shifted somewhat. In our dynamics 
calculations, we shift the triplet PES in energy to best match the closest crossings from the present ab 
initio calculations. We will say more about this shift after the discussion of the singlet PES. 


For the singlet surface, we constructed an analytic representation that included both the ab 
initio results from the MRCI calculations described above and ab initio results from Ref. [42], that 
emphasized the global minimum of the ground electronic state of CO2. To facilitate the determination 
of this new analytical representation, we needed to have as similar quality data from the two sets of 
ab initio calculations. Thus, we computed MRCI energies at 15 geometries near the CO2 minimum that 
were included in the Huang et al. data set. From Ref. [42] we used the ACPF data computed using the 
triple zeta basis. We compared these results at the 15 common geometries with our MRCI results as 
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Figure 6. Scattered plot of 3A’-1A’ spin-orbit coupling matrix element LSz as a 
function of the energy difference between the 3A’ and 1A’ state near the 
crossing seam, within +1 kcal/mol difference in energy. 
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Figure 7. Scattered plot of 13A”-1A’ spin-orbit coupling matrix element LSx as a 
function of the energy difference between the 13A” and 1A’ state near the 
crossing seam, within +1 kcal/mol difference in energy. 
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Figure 8. Scattered plot of 13A”-1A’ spin-orbit coupling matrix element LSy as a 
function of the energy difference between the 13A” and 1A’ state near the 
crossing seam, within +1 kcal/mol difference in energy. 


well as our MRCI+Q results. We found that the root-mean-square (rms) deviation from the average 
difference to be significantly smaller for the MRCI+Q results, so we used them in the new analytic 
representation. The rms deviation from the average difference was 0.47 mEn (approximately 1 
kJ/mol). For this analytic representation, we used the same functional form as for the triplet surfaces, 
with the exception that in the CO+0 asymptote, the CO potential curve was taken to be the minimum 
of CO(*x*) + O({D) and CO(I) + O(3P). The CO(#IT) potential curve was taken from extensive icMRCI 
calculations using the cc-pVTZ basis with two sets of diffuse functions on C and three on O. [43] The 
non-bonding interaction parameters for CO+O and C+Oz2 diabats [41] were fixed at the values found 
for the 3A’ surface, but the modulating/damping function was set to unity. In the bound CO: diabat, an 
eighth order expansion was used. Since the available singlet ab initio data did not cover the full range 
of coordinate space, we had some difficulty ensuring that the final analytic representation behaved 
physically for all geometries of interest yet also gave a satisfactory representation of the ab initio 
data. In the end, we performed a fit with only a fourth order expansion, and this was used to generate 
a set of constraint data points to ensure that nothing untoward happened when using the eighth 
order fit. In the final fit, 8849 icMRCI+Q points, 577 icACPF points, and 55,000 constraint points were 
included. The icMRCI+Q points had uniform weight equal to 577/8849, the constraint points had 
uniform weight 5.77/55000, and the weight for the icACPF point i was 1/max (0.01E), E; — Emin), 
where E,,;, is the energy at the COz2 singlet minimum. As before [41] we checked the percentiles of 
the fitting errors and removed points that seemed to be outliers. We removed 49 icMRCI+Q energies 
and 31 icACPF energies. 


Because of the unequal weights in the least squares fit, it is hard to characterize the quality 
of the fit with a single number. Instead we plot some comparisons. In Figure 9, we plot the fitted 
energy vs. the calculated energy. A perfect fit would lead to points all lying on the line x=y. The zero 
of energy for the plot is the global minimum of the potential. We see small differences between fitted 
and calculated energies over the entire range. For reference, the computed dissociation energy to 
CO(‘=*) + O(‘D) is 255 mE» (675 kJ/mol). Next, in Figure 10 we show the quality of the fit close to the 
minimum. The low lying ACPF energies are fit to within +0.2 mE» while the low lying MRCI+Q 
energies are fit to within about +0.5 mEp. The difference in quality of the fits is perfectly consistent 
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with the 0.47 mE; rms deviation between the two sets of ab initio calculations at the common set of 
geometries. The Harmonic vibration frequencies computed from this surface are @.,=1348 cm-1, 


@p=665 cm-l, and @,s=2362 cm-l, where SS is the symmetric stretch, B is bend, and AS is 


antisymmetric stretch. In comparison, the values from the very accurate PES of Ref.[42] are 1354, 
673, and 2396, respectively. Thus, this is a very satisfactory fit. 


We have to set a reference energy (i.e., a zero of energy) for some COz geometry that is 
accurately represented on all the PESs. Ideally, the zero of energy for the triplet and singlet PESs is 
determined from their separated atom limits, using the experimentally determined O(‘D) - O(?P) 
energy difference. However, this region is not sampled in the ab initio data used for the singlet PES 
determination, so this is not a useful method. Alternatively, we could set the zero of energy by 
making the CO(!x*) + O(‘D) energy relative to the CO(!x*+) + O(3P) equal to the experimental value 
[44] of 15,870 cm. But, as mentioned above, while the analytic representations of the triplet PES 
reproduce well the shape in the vicinity of the crossing seams, there are systematic differences away 
from the crossing seams. We compared the triplet analytical representation to the calculated triplet 
energies at all geometries where the computed difference between the singlet and triplet energies 
was less than 0.1 mE; and we computed the average difference and rms deviation from the average. 
For the 3A’ PES, the average difference and rms deviation for 42 energies was 173.6 and 6.4 mE; 
respectively and for the 13A” surface, the average difference and rms deviation for 22 energies was 
158.1 and 7.9 mE; respectively. There is only one crossing energy satisfying this criterion for the 23A” 
surface: the difference there is 157.6 mE. Because of the 15 mE, difference in shift, it doesn’t make 
sense to couple more than one triplet PES in the dynamics calculations. We will eliminate this 
difficulty in future work by refitting the triplet potentials including the present ab initio results along 
with the original ab initio results used to determine the analytical representation. [41] With this zero 
of energy for the 3A’ PES, the difference in energy between CO(!i+) + O(‘D) and CO(‘=*) + O(3P) is 
58.6 mE» (12,861 cm, i.e. 81% of the experimental value. 


The geometry with the lowest singlet energy, where we compute a crossing, occurs at a CO 
bond length of 1.190 A. This is quite close to the global singlet CO2 minimum with a CO bond distance 
of 1.17 A. In Figure 11 we show a contour plot of the singlet PES for one CO bond length fixed at 1.17 
A. We also show the geometries of the ab initio calculations for one CO bond length being within 0.03 
A of 1.2 A. The deep well is clearly seen and there is a barrier to dissociation of about 16 mE). We 
also show the crossing seam +20 mE». In Figure 12 we show the same set of contours for the 3A’ PES. 
Although there is an attractive region for the 3A’ PES [41], the main reason for the attractive contours 
is because the zero of energy is the same as for the singlet PES, namely the CO(!+) + O(‘D) energy of 
58.6 mEp. 


C. Surface hopping quasi-classical trajectory dynamics 


The surface hopping calculations were carried out using a modification of the trajectory code 
described previously. [45] This is a vectorized variable stepsize trajectory code (VVTC) that has been 
used for many different collision processes. We implemented the fewest switches algorithm of Tully 
[29] with the additional decoherence correction of Granucci and Persico. [46] This is done by adding 
2Ns"f additional coupled equations to the equations of motion for the trajectory, where Ns“"f is the 
number of surfaces coupled. These additional equations govern the time dependence of the real and 
imaginary parts of the amplitude on each surface. Once the zero of potential energy was shifted to be 
at the CO(4=*+) + O($P) minimum, we had no trouble accurately integrating these additional equations. 
This shift was required because although the trajectories themselves do not depend on the zero of 
energy but rather just on the gradient of the energy, the leading contribution to the phase of the 
amplitude is proportional to the current potential energy times time. In our initial calculations, the 
zero of energy was taken to be the fully dissociated separated atom limit, and since the energy of 
CO(‘=*) + OP) minimum is about -190 E, with respect to that zero, we obtained extremely rapid 
phase changes that were hard to track accurately. 
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Our code automatically cuts down the time step when the amplitudes change rapidly, so we 
have no trouble identifying regions where hopping can occur. It is convenient for debugging 
purposes to make trajectories independent of whether they are run individually or in a vectorized 
group. In the unmodified code, this is done by setting up the initial conditions sequentially, so the 
random number string is the same in both cases. The addition of hopping introduces the need for 
random numbers along the trajectory, so the previous strategy is no longer effective. The way we 
solved this problem was to use a different random number generator for hopping than was used to 
set up the initial conditions. Since the hopping random numbers are only used sequentially, the 
random number generator does not need to be as sophisticated as we use for the initial conditions. 
[45] Thus we used the routine ranO from Press et al. [47] with the initial seed generated from the 
random number generator used to generate the initial conditions. 
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Figure 11. The singlet PES with one CO bond length fixed at 
2.201 bohr (1.165 A). One O atom is on the y axis at 1.1005 
and the C atom is on the y axis at -1.1005. The contours 
show the energy when the other O atom is at that position. 
The zero of energy is CO(‘*) + O(‘D), and the gray contour 
line is for this energy. The blue (red) contour lines show 
negative (positive) energies in steps of 50mE;. The green 
contour lines show the crossing seam +20mEn. The blue 
symbols mark the ACPF energy geometries, and the green 
symbols the MRCI+Q geometries. 
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Figure 12. PES for 3A’ at the same fixed CO bond length as in 
Figure 11. 


We carried out trajectories using both the adiabatic electronic basis and the diabatic 
electronic basis. In the adiabatic electronic basis, we diagonalize the electronic Hamilton in the 
diabatic basis whenever the potential energy or its gradients are required. In this case, the coupling 
is determined by the matrix element (K|d/dQ;|N), where K labels the current electronic state and N 


labels a generic electronic state, and Q; is one of the coordinates being following along the trajectory. 
dHnm 


This quantity is given by (K|d/dQ;|N) = — YnmVnx “ao, Vmu/ (Ex — Ey), where n and m label 
diabatic states, V,, is the eigenvector for adiabatic state K,H,,, is the diabatic electronic 
Hamiltonian, and EF, is the adiabatic energy for adiabatic state K. This coupling is quite small except 
in the vicinity of crossing seams, where it can become quite large. Thus, trajectories encountering a 
crossing seam have a very high probability of hopping even if the diabatic coupling is small and the 
physical effect is to follow the diabatic surface. In this case, the surface hopping algorithm corrects 
this sudden hop by hopping back down usually in less than a picosecond. This result, however, is not 
consistent with the philosophy of the fewest switches algorithm. In contrast, if the trajectory is 
carried out in the diabatic basis, the coupling is solely determined by H,,,,, which is taken to be a 
constant, because the magnitude of spin orbit coupling is nearly constant over the range of 
geometries considered. Based on our calculations of the spin orbit matrix elements, we used a value 
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of 100 cm for the present calculations. Thus, the frequency of hopping is much lower than in the 
adiabatic basis but, as the coupling never goes to zero, the electronic state amplitudes show small 
amplitude oscillation even when the collision partners are widely separated. Also, when we correct 
the kinetic energy after a hop [29], we base this correction on (K|d/dQ,|N) even though we follow 
the trajectories in the diabatic basis. Neither the adiabatic or diabatic basis is obviously superior to 
the other, and we choose the diabatic basis for all our calculations reported here. Furthermore, 
because of the reference energy issues described above, we only couple the singlet PES to the 3A’ 
surface. 


We carried out two kinds of calculations: one starting with CO(‘x*) + O(3P), and one starting 
with singlet COz. In the former case, we sampled CO ro-vibrational levels from a Boltzmann 
distribution at 5000K and initial relative translational energies of 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 
and 60 mE, (26.25 to 157.5 kJ/mol). We ran 120,000 trajectories for each initial relative 
translational energy, and integrated each trajectory until either an atom plus a diatom emerged from 
the collision or the time of the trajectory exceeded 10° atomic units of time. This time limit is about 
24 picoseconds. For all collision energies, there were some trajectories that were stopped at this 24 
picosecond limit, but only a handful. Results for a typical initial condition (Erg, = 30 mEn) are shown 
in figure 13, where we plot the number of hops vs. the time delay. The time delay is computed as the 
time of the trajectory minus the time the O would need to reach the center of mass of the CO from the 
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Figure 13. Number of hops vs. time delay for O(?P) + CO(42*) 
relative translational energy 30 mE, (78.8 kJ/mol). 


start of the trajectory in the absence of interatomic forces minus the time the O would require to 
leave the center of mass of the CO to the end of the trajectory in the absence of interatomic forces. 
Due to the repulsive nature of the interaction when O gets close to CO, some time delays are negative. 
For this E’re, most trajectories led to an unstable singlet COz molecule that fell part within about 12 
picoseconds, but 6 were still bound after 24 picoseconds. Thus, the probability of forming a long- 
lived singlet COz molecule is on the order or 10. This results in poor statistics for elucidating the 
lifetimes, thus we carried out the second type of calculations that can be construed as a way to 
directly focus on only the trajectories that formed !COz with any appreciable lifetime. In this 
approach, we study the reverse process, dissociation of energetic COz2 molecules by hopping to the 
lowest-energy triplet PES. 


The specification of the initial conditions for a diatomic molecule are well known, but the 
specification of the initial conditions for a triatomic has posed significant difficulties. The problem 
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relates to the fact that for a diatomic, the vibrational motion is decoupled from the rotational motion, 
only depending on the rotational angular momentum, while the rotational motion is simply driven by 
the vibrational motion. [48] In contrast for a triatomic, all vibrations and rotations are fully coupled. 
The customary generalization to triatomics is the determination of good action-angle variables that 
describe the motion of the triatomic in terms of separable motions, however this is not possible to do 
systematically, for the coupled motion can lead to resonances or chaotic trajectories that so distort 
phase space that the determination of good action-angle variables becomes impossible. In our study 
of the H20 molecule [49], we found that only a few vibrational levels could be determined when 
adding rotation, because chaotic trajectories were found with only small amount of rotational 
angular momentum. 


Thus, in the present work we use a different tactic. We will only specify an initial energy fora 
metastable COz molecule (one with total energy higher than the lowest CO(!=*) + O(3P) asymptote 
and lower than the CO(‘=*) + O(‘D) asymptote), and randomly sample over all possible vibrational 
and rotational motion. We do this by specifying the Cartesian Radau [50] vectors describing the CO2z 
molecule in a laboratory fixed frame of reference and use uniform sampling for each component of 
the vectors within the interval +1.3 A. We evaluate the potential energy at this geometry and if this 
energy is greater than the target energy, we sample again and repeat until an energy less than the 
desired energy was obtained. When position coordinates are obtained that satisfy the constraints, we 
use uniform sampling to obtain the time derivatives of all but one coordinate. The sampling box size 
was set to the maximum allowable time derivative based on the available kinetic energy. The 
remaining time derivative was computed from the square root of the difference between the 
available kinetic energy and the kinetic energy from the other components. If this difference was 
negative, another trial was made until the difference was positive. Then finally we randomly chose 
the sign of the square root term. 


For these calculations, we ran a batch of 100 trajectories for each of the energies used for the 
recombination calculations. We analyzed the results by looking at the times for dissociation, 
computed like the time delay described above. This analysis is only valid if each trajectory is equally 
weighted, which is what we have done in the present work. More careful consideration of the initial 
conditions shows unequal weights will be required but that is beyond the scope of the present work. 
We sort the trajectories in order of the dissociation time, and a typical result is shown in Figure 14 
where we show the results for an initial energy of 30mE;, the same as used for Figure 13. We show 
results for two different maximum integration times: 48 picoseconds and 480 picoseconds. Note that 
the latter integration time is ~1000 times shorter than the mean time between collisions in the 
afterbody flowfield. In both cases, we find trajectories that never dissociated, and of the trajectories 
that did dissociate, we can fit the dissociation percent reasonably well with a half-life of 6 
picoseconds. In this fit, we scaled the lifetime plot so 100% fell at 76%, i.e. 24% of the trajectories 
were assumed unable to dissociate for this analysis. These results seem to be consistent with those 
shown in Figure 13. 


It is not yet clear why 24% of the trajectories were unable to dissociate. There are a few 
different possibilities, but it seems likely that there are at least two different regions of phase-space, 
76% falling into one region and 24% falling into the other, and there is some barrier between them. 
Whether this is a barrier that would manifest itself in accurate quantum calculations, or some sort of 
artifact of classical mechanics, remains a question for further research. Nonetheless, in Table 1 we 
give the half-lives and percentage of dissociation determined from the present work. 
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Figure 14. The dissociation time as a function of percent dissociated for the 
same energy as Figure 13. The results are given for two different maximum 
integration times (48 ps in red and 480 ps in green) are given, along with 
the fitted decay function (dotted line). 


The recombination rate for CO(‘X*) + O(?P) to form singlet COz, was then estimated. To 
properly compute this rate, we need to solve the master equation, and this involves following the 
populations of the various states of COz2 as a function of time, which requires state-to-state rate 
coefficients - this is beyond the scope of the present work. Nonetheless we can make some 
reasonable simplifying approximations: see Refs. [45] and [51] where a comparison is made between 
accurate calculations and the simplified calculations described below. The results were within a 
factor of about 2. To compute the rate coefficient for CO + O + M — COz2 + M, we assume the three- 
body recombination for the process takes place by a two-step mechanism: CO + O = CO2'(i), 
followed by COz2*(i) + M — COz + M. In these reactions, i represents a metastable state, i.e. one with 
energy between the triplet and singlet asymptotes. We approximate all metastable states by the 
levels with the energies for which we have computed lifetimes - see Table I. These computed 
lifetimes satisfy the equilibrium criterion in Ref. [45]. Assuming the first step of this mechanism is in 
equilibrium, the three-body recombination rate is given by k?? = Yin Keq (T)k/_(T), where Kéq is the 
equilibrium constant, and k/ is the de-excitation rate. The sum over i includes all CO2 lying between 
the triplet and singlet dissociation limits, and the sum over n includes all bound states of CO2z. We use, 
as a rough estimate, the value of 28 A2 for the cross section for the de-excitation process summed 
over n. This makes ki (T) equal to 10 times the average velocity. The equilibrium coefficient is given 
by the ratio of the partition functions for CO+O and CO2*(i). The atomic and CO partition functions 
are easily computed, the CO2 partition function is much more complicated. To compute it we need to 
know the density of the CO2*(i) states as well as the total density of states for the COz molecule. We 
will estimate this via a procedure similar to that used to generate the initial conditions of the 
trajectories starting as singlet CO2. This calculation is currently in progress. 


Given the low probability for COz recombination (10-4) coupled with the rather low collision 
frequency predicted for the afterbody flowfield, it seems unlikely that a significant amount of singlet 
COz can be formed by this mechanism. Therefore, other processes should be considered, such as the 
exchange reaction CO + 02 < CO2 + O and stabilization of metastable CO2*(i) by photoemission. 


18 


Downloaded by NASA AMES RESEARCH CENTER on July 7, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-3486 


Table 1. Summary of Classical Trajectory Results for Metastable CO2 


Total energy of CO2 in mE? Half Life in ps % that dissociated 
20 13 63 
25 18 63 
30 6.0 76 
35 4.8 77 
40 2.9 81 
45 2.4 85 
50 2.4 90 
55 2.4 90 
60 1.5 94 


a Energy relative to CO(!x+) + O(°P) 


II. Summary 


In this work, we have focused on the two-body recombination reaction of CO2, in which a 
ground electronic state CO molecule and a triplet state oxygen atom collide and make a transition to 
the singlet potential energy surface of COz. We computed the four relevant potential energy surfaces 
(one singlet and three triplet) and the spin-orbit coupling elements between them for CO2 using ab 
initio quantum chemistry methods. We used the calculated PES data to obtain and calibrate analytical 
expressions for these potentials and we used them for surface hopping classical trajectory 
calculations based on Tully’s fewest switches surface hopping algorithm. The results of the trajectory 
calculations indicate that the probability of a triplet-to-singlet transition is approximately 10. 
Calculation of the recombination rate coefficients, for temperatures characteristic of the afterbody 
flowfields of Mars entry vehicles, are in progress. It is unlikely that this recombination reaction is 
responsible for the larger than expected afterbody heating that has been measured during the Mars 
Science laboratory entry. We are also computing COz2 dissociation rate coefficients (the reverse 
process of recombination) using the same surface hopping method. This work represents the first 
theoretical study of these reactions for COz that incorporates an accurate treatment of the singlet- 
triplet transition. 
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